Ground state of a resonantly interacting Bose gas 
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We show that a two-channel mean-field theory for a Bose gas near a Feshbach resonance allows 
for an analytic computation of the chemical potential, and therefore the universal constant /3, at 
unitarity. To improve on this mean-field theory, which physically neglects condensate depletion, we 
study a variational Jastrow ansatz for the ground-state wave function and use the hypernetted-chain 
approximation to minimize the energy for all positive values of the scattering length. We also show 
that other important physical quantities such as Tan's contact and the condensate fraction can be 
directly obtained from this approach. 



I. INTRODUCTION 

In recent years, ultracold Fermi gases have been ex- 
tensively studied near a Feshbach resonance, both the- 
oretically and experimentally. Only more recently, the 
strongly interacting regime for an atomic Bose gas is be- 
ginning to be explored. The reason for this is that ex- 
periments are troubled by strong inelastic atom losses 
in this case. Nevertheless, a number of groups are now 
starting to carry out experiments in the strongly interact- 
ing regime near a Feshbach resonance [IH3] , and it is ex- 
pected that significant results will soon be obtained. The- 
oretically, these systems are also challenging and some 
attempts toward an accurate description of the strong in- 
teraction effects that play a role here have already been 
made @HHj- In this paper we discuss another approach 
to study the ground state of resonantly interacting Bose 
gases. 

As a first step and to discuss more transparently some 
of the important physics involved, we start this paper in 
Sec. [n] with a mean-field description of an atomic Bose 
gas near a Feshbach resonance. This mean-field theory 
is based on a two-channel description containing both 
atoms and molecules, and has as a main approximation 
the neglect of depletion of the condensate. Using a two- 
channel model gives a finite energy for the Bose gas for all 
values of the scattering length a, also at unitarity, where 
the scattering length diverges. Moreover, near the Fesh- 
bach resonance the theory can be written in a universal 
form, which no longer depends on the specific details of 
the system. In this form it is even possible to find an 
analytic solution for the chemical potential at resonance. 

However, this mean-field theory is not qualitatively re- 
liable for large interaction strengths, since it neglects con- 
densate depletion, which has significant effects on the en- 
ergy. Therefore, we also study in Sec. |III| a variational 
Jastrow ground-state wave function combined with the 
hypernetted-chain approximation for the calculation of 
the ground-state energy. This approach has had great 
success in the strongly coupled helium liquids and we also 
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show that, as desired, it reduces to the Bogoliubov theory 
in the weakly interacting limit. After a somewhat techni- 
cal description in Sees. |III AfjinC] on h ow to implement 
this approach, we show in Sees. Ill D| [IiTF| that it can 
be used to directly compute several important physical 
quantities. For instance, using this approach, the conden- 
sate fraction and the contact can be derived directly from 
the two-particle correlation function. As mentioned, the 
approach is variational. The total energy of the gas can 
be determined from the two-particle correlation function 
and the Jastrow factor, which are related to each other 
via the hypernetted-chain equation. The Jastrow fac- 
tor, which determines the many-body ground-state wave 
function, is ultimately found by minimizing the energy. 
In Sec. IIVI we find that for the small and intermediate 
scattering length regime na^ < 1, where n is the atomic 
density, this approach works very well, and also allows us 
to compute the contact and condensate fraction. How- 
ever, the parametrization of the two-particle correlation 
function that we use here and that is inspired by the liq- 
uid helium literature, does not appear to work properly 
for larger scattering lengths and this remains a topic for 
future work. 



II. MEAN-FIELD THEORY 

In ultracold dilute Bose gases, the interactions are 
usually completely determined by the s-wave scattering 
length a. However, the two-atom scattering problem can 
also contain bound states. In the case of a magnetic 
Feshbach resonance, the energy of these bound states de- 
pends on the externally applied magnetic field B. At 
certain values of this magnetic field, a new bound state 
can cross into the continuum of scattering states. At such 
a point there is a resonance in the scattering length, and 
the interaction appears infinitely strong in the s-wave 
channel. 

In order to describe the many-body physics in such a 
system we start with an effective action for the atom field 
(0a) and molecule field {(j)m) that describes the bound 
state. This action can be derived from first principles 
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where n is the chemical potential, V is the volume, and 
/3 = l/k-aT is the inverse temperature. For our purposes 
we can restrict ourselves to the zero-momentum and zero- 
frequency part. The atoms interact with each other in the 
so-called open channel with a strength Tbg, proportional 
to the background scattering length abg for the atoms, 
which is an experimentally known property of the spe- 
cific Feshbach resonance of interest. The width of the 
resonance is determined by the atom-molecule coupling 
g and is also known experimentally. The molecular en- 
ergy depends on the external magnetic field through the 
self-energy KEi^ and via the magnetic detuning 5{B) ex 
B — Bq from the resonance at the magnetic field Bq. For 
very broad resonances the interaction strength g of the 
atoms with the molecules can, in principle, also depend 
on the magnetic field, but we neglect this feature here. 

Minimizing the action gives rise to the following Gross- 
Pitaevskii equations for the atoms and molecules. 
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where we have introduced the Hartree-Fock self-energy of 
the noncondensed atoms /iShf- Its precise form is given 
below. The introduction of this self-energy is very impor- 
tant. Without the shift of the self-energy of the noncon- 
densed atoms, the molecular condensate would always 
be unstable. Incorporating the Hartree-Fock self-energy 
makes sure that a condensate of molecules does not de- 
cay away immediately. In other words, the Hartree-Fock 
contribution to the self-energy makes sure that there ex- 
ists a (metastable) equilibrium solution of the mean-field 
equations. Note that by elimination of the molecular 
field and considering the two-body limit, it is easy to 
show that the effective T matrix of the atoms obeys the 
standard relation for the scattering length; that is, the 
effective scattering length a{B) is related to the magnetic 
field via 
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For the broad Feshbach resonances of interest to us 
here, the molecular field and therefore the molecular den- 
sity turn out to be very small, and we are allowed to put 
the atom density ria equal to the total density n. As 
a consequence, the two-channel model now reduces to a 
single-channel model. The mean-field theory now reduces 



to solving the following three coupled equations 
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where abg is the background scattering length associated 
with Tbg- The first equation follows from the Gross- 
Pitaevskii equations in Eq. (pi), the second equation is 
the standard form of the molecule self-energy first de- 
rived in this context in Ref. [8] and the third equation is 
the appropriate Hartree-Fock self-energy. 

In this paper we are especially interested in the uni- 
tarity limit, which is the limit a(B) — > 00. The physical 
properties of the atomic Bose gas are in this limit uni- 
versal, which means that these properties do not depend 
on the specific details of the system, such as abg and g. 
This can be seen explicitly from the equation above. In 
the limit that a — > 00 the background scattering length 
flbg is irrelevant. Thus, we are allowed to take the limit 
flbg — >■ 0, while keeping g^ /S{B) constant and still obey- 
ing Eq. ([3]). Furthermore, the experimentally interesting 
case is a broad Feshbach resonance; we therefore take 
the limit g -^ 00 and 5{B) — >■ 00, while keeping the 
scattering length a constant. In order to proceed fur- 
ther we introduce the Fermi momentum fcp and Fermi 
energy ep instead of the density n = kp/Gn'^ and the 
mass m/fi^ — fc|/2eF. We then end up with 
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where we used the following universal relations for a(B) 
and hT.^{E): 



29^ 

5{B) 

5{B) 



8'K{k-pa 



ep 



=k'pa\ — 



E 
2ep' 



The former two equations give the chemical potential and 
the Hartree-Fock self-energy in units of the Fermi energy. 
The two equations in Eq. ([5| can be solved (in prac- 
tice numerically) for any positive value of a. The result 
is shown in Fig.[T] For small a, the relation for the chem- 
ical potential simply reduces to // = 4eFfcFa/37r = nT{a), 
which is the well known Gross-Pitaevskii expression for 
the small a regime. As expected the Hartree-Fock self- 
energy then reduces to /lE^^ = 2nT{a) = 2/d. We can 
also solve the chemical potential explicitly in the unitar- 



ity limit. We then have 
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Here we notice immediately that fiS^^ = 2\/2/x, from 
which we can then easily solve for /i to obtain 
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For fermions, there is a similar relation, which is usually 
written as fi — {1 + /3)ep. In the specific case of fermions, 
the universal constant /3 contains all the interaction ef- 
fects and was found to be /3 ~ —0.58 [5Vll2j. It is custom- 
ary to define a similar /3 for bosons; the above mean-field 
theory gives /? ~ —0.54. This is just below an experimen- 
tal lower bound set at /3 > —0.56 [5]. Other theoretical 
analyses give varying results, namely, /? ~ —0.34 j5] and 
the upper bounds /3 < 1.93 [5 and /3 < -0.20 [5]. It is 
remarkable that the fermionic value of (3 is within these 
bounds: thus, it is not excluded that there is for this 
quantity no difference between fermions and bosons at 
unitarity. This might be anticipated in a one-dimensional 
situation; however, for a three-dimensional gas as consid- 
ered here, this would be an interesting result indeed. 

Another well-known mean-field result for the Bose gas 
energy is obtained from Bogoliubov theory, which is an 
expansion in terms of the diluteness parameter na^. This 
was already derived in the late 1950s in Ref. [T3] and 
reads 
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where e is the energy per particle. The first term is 
the Gross-Pitaevskii (GP) result and the second term 
is known as the Lee-Huang- Yang (LHY) correction, and 
is due to condensate depletion resulting from quantum 
fluctuations. 

The energy in Bogoliubov theory diverges in the limit 
a — >■ oo; this is shown in Fig. [l] as the dashed line, which 
includes the LHY correction. The dotted line is the GP 
result without this correction. However, the energy for 
the mean-field theory of Eq. (|5| stays finite, shown in 
Fig. [l] as the solid line. Of course, the energy of the uni- 
tary Bose system must be finite; thus the mean-field re- 
sult of Eq. (51) describes this behavior correctly. However, 
from Fig. IT can be concluded that quantum fluctuations, 
described by the LHY correction, are not properly incor- 
porated in this theory Thus, although qualitatively cor- 
rect, the mean- field approach described in this section is 
probably not very reliable quantitatively. To improve on 
this we propose a different approach, based on a Jastrow 
wave function and the hypernetted-chain approximation 
which is discussed extensively in the rest of this paper. 



l/kpa 

FIG. 1: (Color online) The chemical potential /i in units of the 
Fermi energy ep as a function of the inverse scattering length 
l/fcpa. The solid line shows the mean-field result from Eq. (Isl. 
The dashed line shows the Bogoliubov result from Eq. l8| 
with LHY correction, while the dotted line is without this 
correction. Our two-channel mean-field approach stays finite 
in the unitarity limit, while the Bogoliubov theory diverges 
in that limit. 



III. JASTROW AND HYPERNETTED-CHAIN 
APPROXIMATION 

In the previous section we have shown that with a 
simple mean-field theory, it is possible to capture the 
qualitative behavior of a Bose gas, where the energy 
stays finite when the scattering length diverges. How- 
ever, this approach is probably not able to predict the 
energy reliably, since it excludes quantum fluctuations. 
We therefore propose an alternative approach in which 
we make a Jastrow ansatz for the wave function and use 
the hypernetted-chain approximation to compute corre- 
lations. This method was applied with great success in 
the field of strongly interacting helium [TH [TS] . 

Since the Jastrow ansatz in combination with the 
hypernetted-chain approximation has been used success- 
fully for some time now, there exists a large amount of 
literature on the subject. However, in the field of ultra- 
cold atom gases, it is not used very often. We believe 
that these methods can be important for this field and 
we therefore briefly summarize the important relations 
and derivations in the sections below. 



A. Jastrow ansatz 

The many-particle wave function can be a very com- 
plicated function of all the particle positions, but in the 
Jastrow approximation it is argued that the dominant 
correlation features are captured by the pair function or 
Jastrow factor f{ri — r2) = f{ri2)- The wave function is 
then. 
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In a homogeneous system this Jastrow factor only de- 
pends on the relative positions. This function goes to 
one on a length scale larger than the interparticle dis- 
tance. 

An important function in this description is the two- 
particle correlation function. It is defined as follows: 

Here n is the density, J dR = / dri . . . dr^ denotes the 
integration over all spatial coordinates, while J dRi2 = 
J dra . . . drjv is the integration over all spatial coordi- 
nates except Ti and r2- 

The energy of a system with a Jastrow wave function 
can be written in terms of the functions / and g. The 
potential energy in terms of the Jastrow wave function is 
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where V{r) is the interparticle potential, which only de- 
pends on the distance between the particles. Using the 
particle-exchange symmetry of the wave function we can 
write this in terms of g{r) as 

(V) = ^n^ J dn dr2 V(ri2).g(ri2) . (11) 

The kinetic energy can be written as 

h^ /di2**(ri,...,rAr)Efvf*(ri,...,rjv) 
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which can, again using the symmetry properties of the 
wave function, be written in terms of /(r) and g{r), 

{T) - -^^ / dridrs gin^^logfin^) . (12) 
2 2m J 

Since we describe a homogeneous system, we can perform 
one more spatial integral, which gives a volume factor V. 
The total energy is now, 

e=^nJdrgir)(vir)^^Vl\ogf{r)^ , (13) 

where e is again the energy per particle. 

In this Jastrow ansatz for the wave function, every- 
thing in the system is determined by the Jastrow factor 
/(r). However, many quantities, like the energy, are di- 
rectly related to the two-body correlation function g{r), 
but unfortunately the relation between g{r ) and /(r) is 
very complicated. This relation [Eq. ( 10 )] contains as 
many integrals as there are particles, which clearly is 
unsolvable analytically. There are many approximation 
schemes to solve it, but many depend on small interac- 
tions, or correlation lengths. However, for bosons near a 
Feshbach resonance, we need an approximation scheme 



where these are large. The hypernetted-chain approxi- 
mation, which is a diagrammatic cluster expansion, has 
proven to work also very well in the strongly interacting 
regime [IS]. After we have established the relation be- 
tween / and g, we can solve for / (or g) by minimizing 
the energy in Eq. (13). 



B. Hypernetted-chain approximation 

With the Jastrow wave function we have a direct re- 
lation between the two-particle correlation function g{r) 
and the Jastrow factor /(r), but this relation contains, in 
the thermodynamic limit, an infinite number of integrals. 
These integrals cannot be solved analytically, but using 
the hypernetted-chain (HNC) approximation we can sys- 
tematically evaluate them. The precise details of HNC 
can be read elsewhere, for example, in Rcf. |16| . but 
in the following we give a short derivation for complete- 
ness' sake to understand better the physics involved and 
introduce some useful notation. 

We start by defining the cluster function h{r) = 
f{r)'^ — 1, which goes to zero quickly for large r, since 
/(r) goes then to one. The two-particle correlation func- 
tion g can then in a natural way be written as a cluster 
expansion in terms of /i, 
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where the normalization constant is irrelevant for the 
discussion and is left out and the relative coordinates 
are defined as Vij = ri — tj. These integrals are now 
written as an infinite sum of clusters of /i, which each 
are a product of any number of /I's. These in turn can 
have different levels of complexity in terms of the inte- 
gration variables. For example, J dr^h{ri'^)h{r2'i) is more 
complicated than J dr^ drj^h{riz)h{r24). The idea behind 
hypernetted chain is to sum over an infinite amount of 
clusters selected by their complexity. When all complex- 
ities are taken into account, we end up with the exact 
result. However, in this paper we stick to the simplest 
set of clusters or diagrams, called nodal diagrams. This 
is referred to as HNC/0. Since this is still a sum of an 
infinite amount of diagrams, the convergence of the ap- 
proximation does not depend on the density or interac- 
tion strength to be small. 

The nodal diagrams are all clusters of h where the 
integral over a series of /I's only connects one h to the 
next. Here 'connect' means that we have an integral 
like J dr3/i(ri3)/i(r32), where ry, 'connects' the two cluster 
functions. We can construct an infinite set of these with 



the following recursion relation: 
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where Af-^'^^ denotes the set of these simple nodal dia- 
grams. 

This set can, in turn, be used to generate an infinite 
amount of composite diagrams, which is simply all pos- 
sible products between all the elements oi J\f^^'^\ This 
we can write as 
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where the numerical factors exactly cancel any double 
counting. This set can be extended even further by 
adding h{rab)exp{J\f'-"^^\rab)), leading to 

/2(r,,)exp [AA("^i)(^a6)] -^f'^'''Hrab) - 1 , 

(17) 

where X^^'^'' is a set of all composite diagrams we can 
make with the set TV'^'''^^ 

A lot more diagrams can be constructed by defining a 
set A/'^°'^^ that obeys Eq. (15) but with h{rab) replaced 
by X'^^'^\rah)- We can proceed naturally, and define a 
X^°''^\rab) that obeys Eq. ([It]) where TV^"'^) is replaced 
by A/'^°'^). We can continue doing this, and in the limit 
where this procedure is followed an infinite number of 
times, we arrive at the following recursion relations: 
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where limi._j.oo 



_^(o,fe) ^ _^(o) and limfe_,oo A'fO''^) = A'(O). 
This latter equation relates the two-particle correlation 
function g to the Jastrow factor /, which is what we 
needed. This selected set of diagrams used to compute g 
is called HNC/0. In order to include more(all) contribut- 
ing diagrams we would have to include also more(all) el- 
ementary diagrams in Eq. pTJ ), in addition to the nodal 
diagrams. However, this HNC/0 approximation contains 
already a lot of important information, as was shown by 
the calculations on strongly interacting helium. 



The relation between / and g in Eq. ( 19 ) can be solved 
for / as 



log/2(r)=logg(r)-AA(0)(r). 



(20) 



The usefulness of this equation follows from the fact 
that the function J\f^^\r) can also be related to the two- 
particle distribution function. To do this we first define 
the structure factor S{k), 
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Note that from the definition of g in Eq. ( 10 1 it follows 
that 5(0) = 0. The integral relations in Eqs. 18 and 
|19| can be written as algebraic equations after a Fourier 
transformation. These equations are then easily solved 
and we get A/'''°'(fc) in terms of S{k), 
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where A/'*-°-'(fc) is the Fourier transform oi Af-^^r). In 
the HNC approximation, the Jastrow factor /(r) is thus 
completely determined in terms of the two-particle cor- 
relation function g(r). 

In this paper, we will vary the function g and then 
calculate the energy. For this we need /, which we can 
calculate using the above equations. Since we need some 
complicated shape for g, it is not possible to analytically 
perform the Fourier transformations. Calculating the en- 
ergy thus involves a few steps. First, when we have a 
g, we calculate the structure factor S{k) by numerically 
Fourier transforming this g. Second, we calculate A/"'-"-' (k) 
and numerically inverse Fourier transform back. Third, 
we calculate / with which we can compute the energy. 



C. Energy minimization 

With the relations that follow from the HNC approx- 
imation, we can write the energy of the system in terms 
of only the two-body correlation function. To get the 
ground-state wave function, we have to minimize this en- 
ergy with respect to this function. We first write down an 
analytic expression for this minimization condition, but 
it turns out to be hard to solve this relation in practice. 
It is much more convenient to numerically minimize the 
energy. 

We have a relation for the energy in terms of / and g in 
Eq. (13 1, and in combination with Eq. (20 1 we can write 



this in terms of g only. Taking the functional derivative 
of the energy with respect to g{r), or more conveniently 
\/g{r), and putting that to zero gives the following dif- 
ferential equation for g 



'2^^' + [^(^^ + """(^^l ^ ^^ = " ' ^^^) 



where a;o(r) is defined as the inverse Fourier transform 
of 



4to 



1)1- 



1 



(24) 



The minimization equation Eq. ( |23[ ) has the form of a 
simple Schrodinger equation for ^/g where wo(r') acts as 
an effective induced potential that takes the presence of 
the entire medium into account. This may seem as a 
simple-to-solve equation, but recall that loq (?■) contains g 
in a very non-linear way. 

Solving this differential equation for g numerically 
turns out to be very hard. Small numerical errors in 
g trigger solutions of the differential equation that are 
not physical, that is, these solutions are not normaliz- 
able. By explicitly varying g to minimize the energy, 
these problems can be circumvented. 



D. Asymptotic behavior 



k. This means that the 1/r^ tail will have a different 
prefactor, but should still be there in the unitarity limit. 
The tail of the trial functions for g, which we use in 
the variational calculation, will therefore be of that form. 
From the prefactor we can determine the speed of sound. 
The Jastrow factor / is completely determined by g 
(and S) and the large-r tail of this function is thus also 
known 
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These limits also tell us something about the large-r be- 
havior of the effective induced potential in Eq. ( 24 1 , 
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This result is consistent with the analytic minimization 



equation in Eq. ( 23 1 , since (when we put V to zero) the 



two limits for both g and uiq exactly solve this differential 
equation. 



Even without minimizing the energy we can say some- 
thing about the shape of the two-body correlation func- 
tion g. Let us first study the case for small scattering 
length kpa. In this regime, the known Bogoliubov dis- 
persion relation can be related to the structure factor 
S{k). This relation follows from the dispersion relation 
from Bijl-Feynman theory, which reads 
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In Bogoliubov theory, the dispersion relation is given by. 
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Here, efc = ffk'^/2m and T{a) = iTrffa/m. When we 
combine Eq. ( 25 1 and Eq. ( 26 ) we get the following for 
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Thus, for small k we have S{k) — hk/2mc with c = 
y/Anan the speed of sound of the medium. 

As was pointed out before, the structure factor is re- 
lated to the two-body correlation function g. Since we 
know the behavior of S{k) for small fc, we can deduce 
the large-r behavior of g. Using the asymptotic Fourier 
transform we get for large r that 
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This result holds only for small a, however, for large a one 
still expects to find a linear dispersion relation for small 



E. Contact 

The small-r behavior of the two-body correlation func- 
tion can be related to what is called the contact, denoted 
by C. The contact was recently derived to be a gen- 
eral feature in strongly interacting Fermi systems by Tan 
[T71 [IH], in a sequence of papers published in 2008. Since 
Tan's derivation is not based on the statistics of the par- 
ticles, it was pointed out by Combescot et al. 19 that 
the relations also hold for Bose statistics and are hence 
applicable to Bose gases as well. The quantity C is part 
of a series of various exact and universal relations which 
therefore also hold for strongly correlated gases. When 
applied to Bose gases. Tan's main theorem, which he calls 
the "adiabatic sweep theorem" , states 
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here e is the energy per particle of the gas. It is striking 
that this is such a simple, exact and universal relation. 
The contact C turns out to be independent of the short- 
range interactions, except for the scattering length a. In 
general, it is a constant which is expected to remain fi- 
nite for all values of the scattering length. Let us con- 
sider the well-known low-density expansion, or Bogoli- 
ubov theory for the ground-state energy Eq. (Is]). When 
we apply Tan's theorem to this energy expression, we find 
an approximation for the contact of a Bose gas for small 
scattering length to first order. 



C = {^-Knaf 



(32) 



It was shown [20] that this equation for C can be derived 
independently of Eq. (31), from which can be seen that 



Tan's relations agree with Bogoliubov theory. 



The relation for the contact in Eq. ( 32 ) does not have 



a finite hmit at a Feshbach resonance when a goes to in- 
finity, since Bogohubov is only valid for small na^ . How- 
ever, the method proposed in this paper does have a finite 
limit. It is in general possible to find C in terms of the two 
particle distribution function g{r), which in the context 
of HNC/0 is of great use. 

For small r, the behavior of the two-particle distri- 
bution function g{r) is dominated by the interaction of 
only two particles, since in a dilute gas, the rest of the 
particles are far away. The function ^/g is therefore pro- 
portional to the two-particle wave function f2{f), which 
is the solution of the two-particle Schodinger equation 
(see Sec. IV A). For small r, but outside the range of the 
interaction, this function behaves as /2(?') ~ 1 — -. The 
two-body distribution function g is thus for small r pro- 
portional to /|, with a proportionality factor we call Z. 
Thus, we have for small r. 



g{r)c^Z\h{rt 



Za^ 



1 



2 
ar 



(33) 



The proportionality constant Z is related to the contact 
through Z = C/ldiT^n'^a^. We thus have for the contact 



n ir 2 2 2 9v) 

C ~ Ion n a 



\.f2{rW 



(34) 



which can be calculated directly from the HNC solutions. 
We calculate the energy as a function of the scattering 
length, and as a result, we are also able to use the original 
expression in Eq. (31 1 to compute C, which we compare 
to the results from Eq. (34). 



In the experiments that a number of groups a trying to 
perform, one of the biggest challenges is the severe losses 
of the particles in the trap. An important consequence of 
the wave function renormalization factor Z, which is re- 
lated to the contact, is that it also affects the three-body 
collision terms, which is what determines the particle loss 
rate of the Bose gas in a trap. This particle-loss is gov- 
erned by the relation 



dn 
di 



L It" 



(35) 



where L determines the loss rate, and the power of n"^ 
reflects that three-body collisions are needed to obey the 
conservation laws. Since the wave function amplitude 
changes with Z < 1 at small distances, the loss rate L 
is multiplied with Z'^ due to many-body effects. Since 
the contact, and hence Z, can change significantly near 
the Feshbach resonance, this will have great effect on the 
losses in experiments. 



F. Condensate density 

An important physical quantity is the condensate frac- 
tion, denoted by ng- This is the density of particles 



which are in the zero-momentum state and form a Bose- 
Einstein condensate. Conversely, there is a density of 
particles which are not in the condensate, due to (quan- 
tum) depletion. This density is typically nonzero even at 
zero temperature, an effect which is solely due to inter- 
actions. In this section, we follow the lines of Ristig et 
al, seeRefs. [HHIS]. 

We first consider the one-body density matrix for the 
system of N bosons, given by 

. ^ , Jd.Ri^*(ri,...,rjv)^(r^,...,rjv) 
"^^"'^=^ /dilr|*(r„...,r.)p ' ^''^ 

where ^ is the wave function for the system. This can 
be written in a convenient notation as 



n{r) = noB '^'^^^ 



(37) 



The one-body density matrix has the well-known prop- 
erties that n(0) — n and n(r — ^ oo) = no, which in terms 
of Q means Q{r — >■ oo) — and uq = ne'^'^^\ Using a 
cluster expansion, similar to HNC/0, this Q{r) can be 
computed. The details of this computation can be found 
in Refs. [21105] , but in the following we give a brief out- 
line of it. 

The most insightful approach to the calculation of Q{r) 
is the method proposed by Feenberg [21]. The ground- 
state Jastrow wave function, written in Eq. ^ is not 
necessarily properly normalized. A trial wave function 
that can be properly normalized, and was proposed by 
Feenberg, is given by 

*W(r„...,r.)=e--(^)"''n/(r,-). (38) 

Here, A is a dimensionless parameter which we need to 
calculate. The one-body density matrix can be written 
in terms of this new trial wave function as 

n(rn') = (39) 

N 

■ dR^\¥''-'\2,...,N)\'l[f{n,)f{n,,). 



ne 



The normalization parameter A can be computed from 
Eq. ( 38 ) by comparing a wave function for N and for iV — 



1 bosons, see Ref. [24J for details. From this comparison, 
A can be calculated in several orders of the previously 
discussed cluster function h{r). Up to second order in 
the cluster function we get 



\^D''^\h]+D^'^\h] + ... = D[h] 



(40) 



where 1)1^1 [h] and 13 '■^l [h] are functionals of h of first and 
second order, given by. 



D^^\h] =n f drh{r), 



(41) 



^'" [h] = y / dr2 dra g{r23 - I)h{n2)h{n3) . (42) 



The expression for the density matrix in Eq. ( |39[ ) can 
be expanded in a similar way. However, instead of the 
cluster function h{r), the radial function ({r) = f{r) — 1 
is used. This function has the similar property that it 
goes (quickly) to zero for large r, and hence we can also 
perform a cluster expansion. Since the normalization was 
computed with four copies of /, and thus second order in 
/i, the density matrix is also computed with four copies 
of / and thus to fourth order in (. 

Again, the precise details of this calculation can be 
found in Ref. [53] , but up to second order in p the result 



IS, 



[21 aPI 

n^ ' — ne exp I 



(2i?W[C] -QW(r))x 

Xexp(27?[2][^]_g[2](^)^ ^ 

which in general can be written as, 

n{r) = nexp(2i?[C] - D\h\ - Q{r)) , 



(43) 



(44) 



with g[il(r) + Q[2l(r) + . . . = Q[r). The function Q{r) 
contains every term that still depends on |ri — rj^j, all of 
which go to zero for r — >■ oo. All constant terms turn out 
to have the same functional form as the normalization 
terms, and can be expressed in the same functional D. 
If we compare Eq. ( 44 ) with Eq. ( 37 ) we notice that it 



has exactly the same form. Thus when we take the limit 
r — >■ oo we get for the condensate density. 



no 



iexp{2D[C]-D[h] 



(45) 



When we insert the expression for D up to second order 
in p, we get 



TiQ = n exp 



dfc 

(2^ 



drC(r)2 



iS{k) - 1) c(fc) 



(46) 



Mkf 



here we have used the Fourier transform of C and h in 
order to get rid of double integrals over r. This expression 
for the condensate density only depends on / and g and 
we are now able to compute this for the minimized results 
below. 



IV. VARIATIONAL SOLUTIONS 

With the HNC approximation for the Jastrow wave 
function, we have an expression for the energy in terms 
of the two-particle distribution function. The ground- 
state g{r) minimizes this energy. In the previous section 
we derived a differential equation for g{r) in Eq. ( 23 1 that 



solves the minimization equation. However, this is a very 
nonlinear equation in g(r), since the effective potential 
a;o(r) in Eq. (24) depends on g(r) in a complicated way. 
This makes solving the differential equation very difficult. 
A variational approach, where we directly vary g{r) to 
find an energy minimum, turns out to work much better. 
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FIG. 2: (Color online) The scattering-length a in units of 
_Rc as a function of the dimensionless interaction strength 
Ce/itV^Rt/m). At certain values of Ce the scattering length 
diverges and the system is at a resonance. 



A. Potential with resonance 

In the unitarity limit, it is expected that the system 
behaves universally; this behavior does therefore not de- 
pend on the exact shape of the interaction potential. This 
gives us the possibility to choose a simple potential that 
is numerically convenient, and also contains a 'Feshbach' 
resonance to go to the unitarity limit. The potential we 
choose is a hard core combined with an attractive 1/r^ 
tail. One of the advantages of this potential is that the 
two-particle problem can be solved exactly. From the 
two-particle solutions, the scattering length can be deter- 
mined, which diverges for certain values of the interaction 
strength of the attractive tail. 

The interaction potential has a hard core with radius 
-Re, and has the following form. 



F(r) = 



-Cs 



r<i?c 
r>Rr 



(47) 



This potential is spherically symmetric and since we are 
looking at dilute and ultracold gases, only the s-wave 
part of the interaction is important. In the spherically 
symmetric case it is convenient to define 



f2{r) = 



u{r) 



where /2(?') is the two-particle wave function. 
Schrodinger equation for u(r) can be written as. 



m dr^ 



^]uir)=0 



This differential equation is solved by 

' \JraC^ \ ( \JmCQ 



u(r) = -sfr 



C\J_ 



2hr^ 



C2J1 



2hr^ 



(48) 
The 

(49) 
(50) 



where Jj.i{r) is the Bessel function of the first kind. 
The hard core is included with the boundary condition 



u{Rc) — and the normalization of the wave function 
demands that f2{r — >■ cx)) = 1. These two relations fix 
the constants Ci and C2 and we obtain 



/2M = 



2hr^ 



(51) 



v^r(f) 



-\mC,/h^)-^/^^rJ,^[^) 



From the two-particle wave function we can determine 
the scattering length a, 



iv{r)J-.{Wf) (C,mV'' 



2r(|) 



ji 






;i2 



(52) 



In Fig. [2] the scattering length is plotted as a function 
of the dimensionless interaction strength CQ/{h?'R'^/m). 
There clearly are resonances at certain values for Cg . We 
have checked that different shapes of the interaction po- 
tential give the same results for the HNC calculation as 
long as the scattering length is the same. This is due to 
the dilute limit in which the interaction is governed by 
the scattering length, and therefore the energy is only a 
function of a and, in that sense, independent of Cg. 



B. Varying the radial distribution function 




k/kv 



FIG. 3: (Color online) (Top) The minimized radial distri- 
bution function g (solid line) and f{rY (dashed line) as a 
function of the radius in units of the interparticle distance 
Ri. (Bottom) The structure factor S{k) (solid line) and the 
dispersion relation in units of ep (dashed line) as a function 
of the momentum k/k-p. 



Within the Jastrow ansatz and the HNC approxima- 
tion, the energy of the system is completely determined 
by the two-particle distribution function g{r). In the 
approach we propose here, we start with an ansatz for 
g{r) that closely enough resembles the expected func- 
tional form, but parametrize enough freedom such that 
we can find, or get very close to, the actual energy mini- 
mum. 

The ansatz for g{r) can be constructed out of three 
parts. The first part is the tail of g, that is, the power 
of r wit h w hich g — 1 approaches zero. As we have seen 
in Eq. (28), g{r) goes to one for large r as 1 — P/^r^'^ , 
where in the weak-coupling limit we also know that 
P4 = h/2Tr'^nmc. The second part of g is the short- 
range regime. Since HNC incorporates the effects of all 
particles onto each other, which in this dilute situation 
is a long-range effect, it has little effect on the short- 
range behavior of the system. It is therefore reasonable 
to assume that for small r, g is proportional to the two- 
body function /2(f)^- The proportionality constant be- 
tween /I and g is related to the contact as discussed in 



Sec. HIE The third part of g, which is left over, is the 
intermediate-range regime. This is where the short- and 
long-range parts are smoothly connected to each other. 
More important for this regime is the normalization con- 
dition of g. This normalization follows directly from the 
definition of g(r) in Eq. (10) and can be written as 



/ 



dr(l-g(r)) = l. 



(53) 



To account for all of this, the middle part needs the most 
variational freedom. The amount of parameters can be 
extended by adding a cosine oscillation to the middle 
part. These oscillations can be expected to be important 
for the strongly interacting regime, where liquid-like shell 
structures may occur, although we have not observed this 
yet. 

Since g{r) is a distribution function, it is always larger 
than zero, and we are therefore able to write it as an 
exponent of another function. Writing it this way has the 
advantage that it cannot accidentally become negative 
when varying the parameters. We split up the ansatz for 
g(r) in a short- (us), a middle- (um), and a long-range 
(ui) part: 



Us = (2 log /2(r) + Pg) X exp (-P^r 



(54) 



«, = -^,-^x(l-cxp(-Pnr-^-)) . 



(56) 



The radial distribution function g(r) is then given by, 



g{r) = exp ( 



ui) . 



(57) 



This parametrization of the radial distribution function 
was common practice in the field of liquid *He, as for 
instance in Ref. H^. 
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FIG. 4: (Color online) The energy as a function of the inverse 
scattering length 1/fcpa calculated with the variational ap- 
proach of HNC/0 (solid line). The (red-shaded) area depicts 
the estimated accuracy of the result, which shows the large 
error for large fcpa. The dash-dotted line shows the Gross- 
Pitaevskii energy while the dashed line also includes the LHY 
correction in Eq. ([8|. This shows that the HNC calculation 
correctly includes this term. 



Ijk-pa 

FIG. 5: (Color online) The condensate fraction as a func- 
tion of the inverse scattering length l/fcpa calculated with the 
variational approach of HNC/O (solid line). The (red) dashed 
line shows the condensate fraction for Bogoliubov theory. The 
condensate fraction for HNC/O is comparable to the Bogoli- 
ubov result for small fcpa, but for larger scattering lengths it 
is significantly smaller. 



C. Results 

In the previous sections we showed how a Jastrow wave 
function, together with the HNC approximation, can be 
used to compute several properties of a Bose gas toward a 
Feshbach resonance. In this section we will show the first 
promising results for small and intermediate scattering 
lengths. 

To find the energy minimum we vary the parameters in 
the distribution function g in Eq. ( 57 ) and use the HNC 
approximation to compute the energy. For the varia- 
tion we use a gradient algorithm which converges slowly 
toward the energy minimum. For small k-pa this goes 
relatively fast and easy, but with increasing k-pa it be- 
comes increasingly difficult. We therefore increased the 
scattering length step by step, and used the resulting pa- 
rameter values of one minimization as a starting point 
for the next. 

In Fig. [3] we show in the top panel the result of a two- 
particle distribution function for which the energy is min- 
imized. Notice the wiggle near r = 0, which shows that 
we are actually dealing with a meta-stable many-body 
solution of the used potential, that acts as the ground 
state in the HNC approximation. We also show /(r)^. 
In the bottom panel the structure factor is shown (solid 
line), which is zero for fc = 0, as it is supposed to be. 
It also starts linearly, and as a result, the dispersion re- 
lation (dashed line) also starts linearly for small fc, but 
becomes of the usual quadratic shape for larger fc. 

This method works excellently for small and also for in- 
termediate scattering lengths. This can be seen in Fig. |4] 
where the solid line shows the energy as a function of k-pa. 
For small kpa {kpa < 0.2), the energy agrees with the 
mean-field result in Bogoliubov theory see Eq. (l8| (dash- 
dotted line). When we increase kpa (0.2 < kpa < 0.5) the 



energy also includes the LHY correction (dashed line). 
When kpa becomes even larger, it becomes increasingly 
harder to find a reliable energy minimum. To indicate 
this, we have estimated the accuracy of the energy. 

Now that we have the two-particle distribution func- 
tion as a function of the scattering length that minimizes 
the energy, we can calculate several other physical quan- 
tities. One such quantity is the condensate fraction. In 
Eq. ( [46| we showed how this condensate fraction can 
be calculated given the radial distribution function. In 
Fig. [5] the condensate fraction is plotted as a function of 
the inverse scattering length l/fcpa. The blue solid line is 
the result from HNC/O and the red dashed line is the re- 
sult from Bogoliubov theory. The result from HNC/O is 
comparable to the Bogoliubov result for small scattering 
lengths, but for larger values of kpa the depletion in the 
HNC/O case is significantly higher than for Bogoliubov 
theory. 

The contact, which was discussed in Sec. |III E| is also 
an important physical quantity. We showed two ways to 
extract the contact from the HNC/O results: one directly 
from the two-particle distribution function Eq. (34); the 
other as a derivative of the energy Eq. (31 1. Since the 



energy follows the Bogoliubov energy, the contact com- 
puted as the derivative of the energy with respect to 1/a, 
is roughly the same as the Bogoliubov contact in Eq. ( 32 1 . 



However, since convergence is not properly reached for 
some of the large-fcpa points in Fig. |4j the contact can- 
not be computed there either. Furthermore, since one 
expects the energy to be finite at unitarity, the contact 
should also be finite in that regime. 

The second method, where we directly read off the 
contact from g, might indicate already that the contact 
becomes smaller than the Bogoliubov result, which is 
shown in Fig. [61 Also, the slow convergence of the vari- 
ational process prevents us from computing the contact 
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1/kpa 

FIG. 6: (Color online) The contact as a function of the in- 
verse scattering length l/kpa calculated with the variational 
approach of HNC/0 (solid line). The (red) dashed line shows 
the contact for the Bogoliubov theory in Eq. (321. 



up to the unitarity limit, but the results for intermediate 
kpa show a decrease in C. The fact that the contact is 
smaller also indicates that the wave function renormaliza- 
tion constant Z is smaller than one. This would indicate 
that the three-body particle decay rate is suppressed by 
many-body effects in the unitarity limit. 



V. CONCLUSION 

The unitary regime for bosons is still not completely 
understood, both experimentally and theoretically. In 
this paper we believe to have shown that the use of a Jas- 
trow ansatz with the HNC approximation gives promis- 
ing results that will help with the understanding. 

In the first section we put forward a very elegant mean- 
field theory which describes a Bose gas near a Feshbach 
resonance. This theory shows the universal nature at a 
Feshbach resonance and can be used to calculate numer- 
ically the chemical potential as a function of the scatter- 
ing length. Moreover, at unitarity this theory gives an 
analytic result for the chemical potential. 

This mean-field theory is probably, for large interac- 
tion strengths, quantitatively not reliable, since it lacks 
the important contributions of quantum fluctuations. 
We therefore propose to use a Jastrow ansatz together 
with the HNC approximation to further investigate the 
strongly interacting Bose gas. We have shown how to 



set up such an approach. From the two-particle distribu- 
tion function g, which is computed with the hypernetted- 
chain relation, several important physical quantities can 
be derived. Not only the energy, but also the condensate 
fraction and the contact can be computed directly from 

9- 

The system of relations for / and g can be solved us- 
ing a variational approach. The two-particle distribution 
function is varied, until the energy is minimized. This 
gives promising results for small and intermediate val- 
ues of the scattering length kpa. For larger fcpa, the 
chosen parametrization of g does not converge in a sta- 
ble manner. It is yet unclear whether this is purely a 
numerical problem or if there are real physical instabili- 
ties involved. Further work is needed to fully understand 
this issue. However, for the regime where convergence is 
reached, we were able to derive the energy, the conden- 
sate fraction, and the contact. 

The ultimate goal would obviously be to find the en- 
ergy exactly at unitarity. To give a first estimate, we 
show in Fig. |4] with the dashed blue and red lines an ex- 
trapolation of the energy. From the mean-field result in 
Sec.lnj we notice that the energy leaves linearly in l/k-pa 
from unitarity. This also corresponds to a constant con- 
tact at unitarity, which is related to the slope of the en- 
ergy. When we assume this behavior to be correct, we 
find an energy e ~ O.Sef. At unitarity the chemical po- 
tential is related to the energy via fi — 5e/3, which results 
in an estimate of the universal number /3 ~ —0.2. This 
is higher than the /3 = —0.54 found using a mean- field 
theory; however, since the convergence for the higher kpa 
energy could not be reached completely, we expect this /3 
to be an upper bound. This is in agreement with the cur- 
rent experiments and calculations [3HB]. For the contact 
we find at unitarity C ~ 10.3 n^''^, which is remarkably 
close to the contact for unitary fermions, C ~ 11 ri"'/^ 
\TE\ . This poses again the interesting question whether 
the universal behavior of fermions and bosons is identical 
at unitarity. 
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